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A simple method is presented for removing the amplitude, frequency and phase modulations from 
the Laser Interferometer Space Antenna (LISA) data stream for sources at any sky location. When 
combined with an excess power trigger or the fast chirp transform, the total demodulation procedure 
allows the majority of LISA sources to be identified without recourse to matched filtering. 
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A revolution is underway that will transform the field 
of astronomy. An array of ground and space based grav- 
itational wave detectors are in various stages of design, 
construction and operation. The first space based detec- 
tor, the Laser Interferometer Space Antenna 0, is being 
developed for a launch in 2011. These gravitational wave 
detectors will open a new window on the Universe that 
complements traditional forms of astronomy . 

However, before we can talk about gravitational wave 
astronomy, a method must be found to turn a gravita- 
tional wave detector into a gravitational wave telescope 
- that is, a device that is capable of locating individual 
sources on the sky and determining their physical char- 
acteristics. In traditional electromagnetic astronomy the 
conversion of a detector - such as the human eye or a 
photographic plate - into an observational instrument is 
achieved by building a device that collects and focuses 
the incident radiation onto the detector. Such an ap- 
proach is impractical for gravitational radiation as the 
weak coupling of the radiation to matter prevents the 
construction of a physical telescope. Instead, a gravita- 
tional wave detector is transformed into a gravitational 
wave telescope by data analysis algorithms. 

In principle we would like to use matched filtering of 
the data against suitable waveform templates as this de- 
livers the highest fidelity results, but in practice the com- 
putational cost of matched filtering can be prohibitive. 
The solution is to adopt a hierarchical approach where 
a fast, but sub-optimal, method is used to find a good 
initial solution, which can then be refined using matched 
filtering. The excess power statistic used in ground 
based gravitational wave astronomy is a good example 
of this approach. However, we cannot simply take the 
ground based data analysis algorithms and apply them 
to the space based systems as the high frequency, ground 
based detectors primarily look for short lived, burst like 
signals, whereas the low frequency space based detectors 
primarily look for long lived, continuous sources. One 
of the key differences between ground and space based 
gravitational wave detection is ratio of the observation 
time to the orbital period of the detector. The ratio is 
small for most high frequency sources, and of order one 
for most low frequency sources. This means that LISA 
data analysis algorithms have to deal with orbital signal 
modulation and multiple overlapping sources. 

Here we describe a simple method for converting the 
output of a space based gravitational wave detector, such 
as LISA, into useful astronomical information about the 



location and physical characteristics a gravitational wave 
source. The method works by removing the amplitude, 
phase and frequency modulation imparted by the orbital 
motion of the detector. For monochromatic sources, the 
demodulation procedure focuses the signal into a narrow 
frequency band, which creates a spike in the demodulated 
power spectrum. By searching for sky locations that yield 
the largest spikes it is possible to locate the brightest 
sources on the sky. These can then be subtracted and 
the procedure iterated until all resolvable sources have 
been found. For chirping sources the demodulation can 
be followed by a fast chirp transform Q , which effectively 
re-concentrates the power that is spread by the frequency 
evolution of the source. Once again the brightest sources 
on the sky can be identified and iteratively subtracted to 
expose the weaker sources. 

The demodulation procedure is remarkably simple. It 
utilizes the fact that space based detectors such as LISA 
are able to return several interferometer outputs A 
general gravitational wave coming from a source in the n 
direction can be expressed in barycentric coordinates as 

h M „(t, h) = h+(t, n)e+„(n) + h x (t, n)e* v (n) (1) 

where h+(t,h) and h x {t,n) describe the two polariza- 
tion states with basis tensors e+,(n) and e*(n). The 
polarization states h + , h x are related to the principal 
polarization states h+, h x by the polarization angle ip: 



h + = h+ cos 2ip + h x sin 2ip 
h x = h x cos 2ip — bK_ sin 2i(j . 



(2) 



In the idealized limit of a single low frequency source and 
a noise-free detector, the interferometer outputs s(t) can 
be expressed as 

s(t) = h + (t,n)D + (t,h) + h x (t,h)D x (t,h) (3) 

where D + (t, h) and D x (t, n) are the sky-location depen- 
dent interferometer response functions |(j and h+(t, h) 
and h x (t,n) are the Doppler shifted gravitational wave 
strains 



R 

h(t,n) = hit H sin#cos(27r/ m i - 

c 



(4) 



Here n — » (6, <fi) denotes the location of the source in 
ecliptic coordinates, R is the Earth-Sun distance and 
/m = 1/year is the modulation frequency. A simple, 
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but impractical, approach for demodulating the signal is 
to form the combinations 



h x (t,n) = 



Sl (t)D+(t,n) - s u (t)D+(t,n) 
D*(t,n)D+(t,h) - D*(t,h)D+(t,n) 



and 



h+(t,n) 



Sl (t)D*(t,n) - s u (t)D? (t,n) 
D( (t, n)D+ (t, n) - D* (t, n)D+(t, h) 



(5) 



(6) 



from interferometer outputs I and II, then perform the 
coordinate transformation 



R 

t = t' sin 8 cos(2ir f m t' 

c 



(7) 



to remove the leading order Doppler modulation Q, 
thereby recovering the demodulated strains h+{t') and 
h x it'). The main problem with this approach is that the 
denominators in Eqns. (0 and © vanish at certain times 
of the year for sources within 60 degrees of the ecliptic 
plane. When detector noise or other sources are present 
the inversion becomes singular. Although this simple de- 
modulation procedure does not work in the time domain, 
it can be made to work in the frequency domain. 

The key advantage of moving to the frequency domain 
is that the singularities in the temporal inversion do not 
correspond to singularities in the frequency space inver- 
sion. By Fourier transforming the data over a finite time 
interval, the demodulation procedure can be re-cast as 
a linear algebra problem that can be solved using the 
robust technique of singular value decomposition. 

To simply our expressions we assume that the observa- 
tion time T is equal to one year so that the frequency res- 
olution A/ = 1/T is equal to the modulation frequency 
f m (the general expressions are not hard to derive, but 
they are more complicated). The Fourier space equiva- 
lents of J3J and (0J become 



s k = A+ (n) h+(n) + A* (n) hf (ft), 



(8) 



and 



h+ = B ln {n) h+{h), hf = B ln (h) h*{n), (9) 
respectively. The matrices A and B are given by 



At l (n)=D+_ l (n), A*(n)=D*_ l (n) (10) 



and 



Bin(n) = J, 



i(;-n)(7r/2-</>) 



(11) 



Here J q (x) is a Bessel function of the first kind of order 
q and we have used the Einstein summation convention 
for repeated indices. It is important to note that the 
matrices A and B are both band diagonal as D~^{h) w 
w D?(h) for \j\ > 4 and J q (x) w for \q\ > x. The 
next step is to interleave the Fourier coefficients of the 



two polarizations h + (t) and h x (t) into a single column 
vector h such that hi n = h+ and /i2n+i = h*. Similarly, 
the Fourier coefficients of the two interferometer outputs 
s\(t) and sn(i) are interleaved into a single column vector 
s. The detector response can then be written as 



s = M(n)h(n) . 



(121 



where the modulation matrix M(n) is the interleaved 
product of A + (n), A x (n) and B(n): 

M 2 k2i = k jBji, M 2 fe2i+i = A* k jBji, 



If the interferometer outputs are sampled N times in a 
year, M(n) will be a 2N x 2N dimensional, complex, 
band diagonal matrix. Multiplying s by the inverse of the 
modulation matrix demodulates any sources located in 
the n direction. In practice, s will contain contributions 
from instrument noise n, and multiple sources \i l (n l ): 



s = ^]VP(n l )h l (n J ) + n. 



(14) 



We can estimate the contribution to the signal s from a 
source located in the n direction by solving for the vector 
h cff (n) that satisfies the equation 



s = M(n)h eff (n 



(15) 



The best fit solution is found using a singular value de- 
composition. Suppose that source j happens to lie in the 
n direction. Then 



,eff/ 



■E 



M~ 1 (ft)M i (n% i (n i ) + M- 1 (fi)n. 

(16) 

The key to the method is that sources at other sky lo- 
cations are not demodulated, and the noise is merely 
reshuffled, not amplified. Thus, we can use the solu- 
tion vector h (n) as an estimate for the source vector 
h 3 (n) . Another key point is that the modulation matrix 
is band diagonal, which allows us to perform the demod- 
ulation across limited bandwidths using small modula- 
tion sub-matrices. If this were not the case we would be 
confronted with the problem of inverting a ~ 10 8 x 10 8 
matrix. Moreover, a source only has to compete with 
other sources that overlap with it in frequency space. 

The demodulation procedure works for any type of 
gravitational wave source, but it is especially useful when 
the source is a mildly eccentric, nearly Newtonian bi- 
nary as the demodulated signal is then approximately 
monochromatic. We can search for these sources by look- 
ing for spikes in the the power spectrum of h cff (n) for 
different sky locations n. A similar approach was em- 
ployed in Ref. using Doppler demodulation of a single 
interferometer output. Since the vast majority of LISA 
sources are non-chirping, circular Newtonian binaries, 
the combined demodulation, power spike search goes a 
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long way toward solving the LISA data analysis prob- 
lem. After establishing the sky location and frequency 
of a source, the next step is to refine the frequency mea- 
surement and extract the Fourier amplitudes of the two 
polarizations h + and h x . Unless the source frequency / 
is an integer multiple of the sample frequency A/ = f m , 
the power will be shared by several discrete Fourier co- 
efficients: 

h+ ~ h + sinc(irx n )e t7TXn where x n = f/f m - n, (17) 

and similarly for h* . Now suppose that the integration 
period T is changed slightly to T' = T(l+e) where e«l. 
The Fourier coefficients of h(t) evaluated over the interval 
T' are related to those evaluated over the interval T by 



h' n = ^ hjsinc(irxj n )e l 



where 



(18) 



(19) 



If the source frequency / is an integer multiple of the 
new sample frequency Af = 1/T', all the power will be 
concentrated in a single Fourier mode h' k . Thus, the fre- 
quency spreading can be removed by finding the value of 
e that maximizes the spike in the stretched power spectra 
of h + (t) and h x (t). This procedure yields an improved 
estimate of the source frequency, and a direct estimate 
of the Fourier amplitudes of the gravitational wave: 



~h+ ~ h+ k e 



i8 



and h x ~ h y 



(20) 



The quantity 6 = 2nf(T' — T) accounts for the phase 
rotation imparted by the change in integration period. 
Once the Fourier amplitudes of the two polarizations 
have been estimated they can be used to calculate the 
amplitude A, inclination t, polarization angle ip and or- 
bital phase (po of the binary using the relations 



, if (K + hx 
ip — - arctan — r— 

4 v Wx-w. 



if (h'-h* 
<Po = t: arctan ^— ^— 

2 v Uf+^ 7 x. 



A = 



arccos 



arctan 



arctan 



A\-A\ 




a + + Ja\-a\ 



(23) 



(24) 
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FIG. 1: The RMS strain spectral density in LISA channels si 
(upper panel) and sn (lower panel) in units of 1CF 22 . 



where 



A x = 



4(cos 2tp hf - sin 2ip h%) 



cos ipo 



4(sin2V) /i^ + cos2-0 h%) 



sin ipo 



(25) 



and ft« = fft(h + ), h{ = $S(h+), h* = K(/i x ), and h{ = 
3?(/i x ). 

As a typical example of a LISA source, consider a 
White Dwarf binary with orbital period 1000 seconds, 
comprised of two O.5M0 White Dwarfs, located 2 Kpc 
from Earth in the direction of the galactic center (6, <j>) — 
(95.54°, 266.83°). Setting V = 111-2°, t = 60° and 
(fa = 20.1°, and using the standard LISA noise curve, 
we arrive at the signal strengths shown in Fig. 1. The 
raw signal to noise SNR r , and the optimal signal to noise 
that can be achieved by matched filtering, SNR m , are de- 
fined: 



SNR r = 
SNR m = 



ibw s h df 



1/2 



(i 2TS h df \ 

VAw &n J 



1/2 



where is the power spectral density of the signal and 
S n is the power spectral density of the noise, and the inte- 
gration is across the bandwidth of the signal. Using these 
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FIG. 2: The RMS strain spectral density of the demodu- 
lated plus h+ (upper panel) and cross h x polarizations (lower 
panel) in units of 10 -22 . 

definitions we find that SNR* = 5.4, SNR* 1 = 5.2, and 
the optimal signal to noise is SNR m = 49.9. Applying 
the demodulation procedure yields the two gravitational 
wave polarizations shown in Fig. 2. The source param- 
eters are recovered to an accuracy of Af = 0.008/ m , 
AO = 1.4°, A(j) = 0.2°, AA/A = 0.054, At = 0.04°, 
Aip = 16.4° and Aipo = 8.6°. These uncertainties are 



comparable to those found using matched filtering, so 
the demodulation method is not only extremely fast, but 
also close to optimal. 

The method continues to perform well when the SNR 
ratios are smaller, or when there are multiple sources in 
the same frequency band. Even with multiple sources 
the procedure delivers a good initial fit for the parame- 
ters describing each source. By combining the demodula- 
tion procedure with the gCLEAN algorithm described in 
Ref. [a|, it will be possible to extract thousands of indi- 
vidual sources from the galactic background at a modest 
computational cost. 

The simple power spike search will not work for rapidly 
evolving systems, such as supermassive black hole bina- 
ries, as even the demodulated power spectrum will be 
spread over thousands of frequency bins. If one thinks of 
a Fourier transform as a matched filter using sines and 
cosines, we see that we get a good match for monochro- 
matic sources (the power spike), but a poor match for 
chirping sources. The solution is to employ a different 
kind of transform that uses filter functions that evolve in 
amplitude and frequency - the Fast Chirp Transform Q . 
Just as a Fast Fourier Transform (FFT) allows us to de- 
tect signals with a constant frequency, the Fast Chirp 
Transform (FCT) allows us to detect signals with vari- 
able frequency. By applying a FCT rather than a FFT 
to the demodulated signal, we can search for the sky lo- 
cations that harbour the brightest chirping sources and 
extract their physical characteristics. 

Much remains to be done in the field of LISA data anal- 
ysis, but fast non-template based methods such as the 
total demodulation procedure and the Fast Chirp Trans- 
form will likely play a key role in future developments. 

The work reported here is supported by the NASA EP- 
SCoR program through Cooperative Agreement NCC5- 
579. 
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